
%%% The Ambiguity of Fishing for Fun, October 2022
%%% This file estimates the parameters for the linear (clogit), CARA, and
%%% alpha-maxmin utility models

%Importing the dataset generated by data_prep.do
results=importdata('ReadyDataFluke.txt'); 
[r1,c1]=find(strcmp(results.textdata,'opt_out'));
IJ=results.data(:,c1);
[r2,c2]=find(strcmp(results.textdata,'chosen'));
I=results.data(:,c2);
[r3,c3]=find(strcmp(results.textdata,'ID'));
ID=results.data(:,c3);
[r4,c4]=find(strcmp(results.textdata,'Priceline'));
C=results.data(:,c4);
[r5,c5]=find(strcmp(results.textdata,'Age'));
[r6,c6]=find(strcmp(results.textdata,'Male'));
Y=results.data(:,c5:c6);
[r7,c7]=find(strcmp(results.textdata,'ReleaseOther'));
[r8,c8]=find(strcmp(results.textdata,'KeepOther'));
Z=results.data(:,c7:c8);
[r9,c9]=find(strcmp(results.textdata,'SFrelease_0'));
[r10,c10]=find(strcmp(results.textdata,'SFrelease_25'));
Q=results.data(:,c9:c10);
[r11,c11]=find(strcmp(results.textdata,'SFkeep_0'));
[r12,c12]=find(strcmp(results.textdata,'SFkeep_25'));
G=results.data(:,c11:c12);
[r13,c13]=find(strcmp(results.textdata,'Pcatch_0'));
[r14,c14]=find(strcmp(results.textdata,'Pcatch_25'));
P=results.data(:,c13:c14);
[r15,c15]=find(strcmp(results.textdata,'SFcaughtmax'));
R1=results.data(:,c15);
[r16,c16]=find(strcmp(results.textdata,'SFcaughtmin'));
R2=results.data(:,c16);
[r17,c17]=find(strcmp(results.textdata,'SFkeptmax'));
R3=results.data(:,c17);
[r18,c18]=find(strcmp(results.textdata,'SFkeptmin'));
R4=results.data(:,c18);
R=horzcat(R1,R2,R3,R4)./100;
K=size(P,2);

rng(4560); 
%Initializing parameters
a=-10;
b=10;
X0L=a+(b-a)*rand(1,13);

%1. Linear Model
myObjL=@(parameters)likelihoodLinear(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters); 
optionsL=optimset('MaxFunEvals',25000,'TolX',1e-14,'MaxIter',25000, 'Display','on');
[xL,fvalL,exitflagL,outputL]=fminsearch(myObjL,X0L,optionsL);
[XL,fvalueL,exiflagLunc, outputLunc, gradL, hessianL]=fminunc(myObjL,xL,optionsL);
normgradL=norm(gradL); 
H_L=-hessianL;
IM_L=-inv(H_L); 
SE_L=diag(sqrt(IM_L))'; 
StatisticL=XL./SE_L; 
conditioningL=cond(H_L); %Checking the conditioning condition ("Verifying the Solution from a Nonlinear Solver: A Case Study", B.D. McCullough and H.D. Vinod, AER, Vol.93, No. 3., 2003)
if conditioningL<(6.7*(10^7))
condition_L='passed';
else
condition_L='failed';
end
%Computing McFadden's Pseudo-R2, AIC and BIC
LL0=21609.7; 
Pseudo_R2_linear=1-(fvalueL/LL0);
LL_l=fvalueL;
aic_linear=2*length(XL)+2*LL_l;
bic_linear=length(XL)*log(size(results.data,1))+2*LL_l;

% 2. Alpha-maxmin Utility Model
%Initializing parameters
X0=[XL(1),XL(2),XL(3),XL(4),rand,XL(5),XL(6),XL(7),XL(8),XL(9),XL(10),XL(11),XL(12), XL(13),rand];

myObj_ambiguity=@(parameters)likelihoodGrad_ambiguity(R,Z,Y,C,I,IJ,ID,parameters);
options=optimset('MaxFunEvals',30000,'TolX',1e-20,'MaxIter',30000, 'Display','iter','GradObj','on');
[X1,fval1,exitflag1,output1]=fminsearch(myObj_ambiguity,X0,options);
[X2,fval2,exiflag2, output2, grad_ambiguity, hessian_ambiguity]=fminunc(myObj_ambiguity,X1,options);
H_ambiguity=-hessian_ambiguity; 
IM_ambiguity=-inv(H_ambiguity); 
SE_ambiguity=diag(sqrt(IM_ambiguity))'; 
Statistic_ambiguity=X2./SE_ambiguity; 
normgrad_ambiguity=norm(grad_ambiguity); 
H_eigenvalues_ambiguity=eig(H_ambiguity); 
conditioning_ambiguity=cond(H_ambiguity);   
if conditioning_ambiguity<(6.7*(10^7))
condition_ambiguity='passed';
else
condition_ambiguity='failed';
end
%Computing McFadden's Pseudo-R2, AIC and BIC
Pseudo_R2_ambiguity=1-(fval2/LL0);
LL_ambiguity=fval2;
aic_ambiguity=2*length(X2)+2*LL_ambiguity;
bic_ambiguity=length(X2)*log(size(results.data,1))+2*LL_ambiguity;

%3 CARA Utility Model
%Initializing parameters
X0=[XL(1),XL(2),XL(3),XL(4),rand,XL(5),XL(6),XL(7),XL(8),XL(9),XL(10),XL(11),XL(12), XL(13)];

myObj_cara=@(parameters)likelihoodGrad_cara(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters);
options=optimset('MaxFunEvals',30000,'TolX',1e-20,'MaxIter',30000, 'Display','iter','GradObj','on');
[X3,fval3,exitflag3,output3]=fminsearch(myObj_cara,X0,options);
[X4,fval4,exiflag4, output4, grad_cara, hessian_cara]=fminunc(myObj_cara,X3,options);
H_cara=-hessian_cara; 
IM_cara=-inv(H_cara);
SE_cara=diag(sqrt(IM_cara))'; 
Statistic_cara=X4./SE_cara; 
normgrad_cara=norm(grad_cara); 
H_eigenvalues_cara=eig(H_cara);
conditioning_cara=cond(H_cara); 
if conditioning_cara<(6.7*(10^7))
condition_cara='passed';
else
condition_cara='failed';
end
%Computing McFadden's Pseudo-R2, AIC and BIC
Pseudo_R2_cara=1-(fval4/LL0);
LL_cara=fval4;
aic_cara=2*length(X4)+2*LL_cara;
bic_cara=length(X4)*log(size(results.data,1))+2*LL_cara;

%Output Tables
dataL=vertcat(XL,SE_L,StatisticL);
data_ambiguity=vertcat(X2,SE_ambiguity,Statistic_ambiguity);
data_cara=vertcat(X4,SE_cara,Statistic_cara);
Table1=horzcat(vertcat(horzcat("SFrelease", "SFkeep", "Otherrelease", "Otherkeep","Risk Aversion", "Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate","Avidity","Male","Ambiguity Aversion"),data_ambiguity),horzcat(vertcat("LogL","Pseudo R2","AIC","BIC"),vertcat(-LL_ambiguity,Pseudo_R2_ambiguity,aic_ambiguity,bic_ambiguity)));
TableA5=horzcat(vertcat(horzcat("SFrelease", "SFkeep", "Otherrelease", "Otherkeep","Risk Aversion", "Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate","Avidity","Male"),data_cara),horzcat(vertcat("LogL","Pseudo R2","AIC","BIC"),vertcat(-LL_cara,Pseudo_R2_cara,aic_cara,bic_cara)));
TableA6=horzcat(vertcat(horzcat("SFrelease", "SFkeep", "Otherrelease", "Otherkeep", "Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate","Avidity","Male"),dataL),horzcat(vertcat("LogL","Pseudo R2", "AIC","BIC"),vertcat(-LL_l,Pseudo_R2_linear,aic_linear,bic_linear)));

%Saving output that will be used in the simulations
save('paramlinear.txt','XL','-ascii'); save('parameters_ambiguity.txt','X2','-ascii'); save('parameters_cara.txt','X4','-ascii');
save('sigmaL.txt', 'IM_L', '-ascii');  save('sigma_ambiguity.txt', 'IM_ambiguity', '-ascii');  save('sigma_cara.txt', 'IM_cara', '-ascii'); 
